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Synopsis 


The present thesis is a study of some algorithmic aspects of /^-skeletons. The 
concept of a /3-skeleton was introduced by Kirkpatrick and Radke [KR85] as a 
way of capturing the “internal shape” of a planar set of points. Indeed, we get 
an entire spectrum of shapes by varying the parameter /?. For a fixed value of /3, 
the “internal shape” is a geometric graph, obtained by joining each pair of points 
whose /^-neighborhood is empty. 

For /3 > 1, the /3-neighborhood of a pair of points pi,pj is the interior of the 
intersection of two circles of radius P\pipj\/2, centered at the points (1 — p/2)pi + 
{P/2)pj and (/3/2)pi -I- (1 — P/2)pj, respectively. For /3 G [0,1], it is the interior 
of the intersection of two circles of radius \piPj\/{ 2 P), passing through p,- and pj. 
We denote this neighborhood by lune0(pi,pj). The parameter /3 characterizes the 
neighborliness of a pair of points. 

Let P = {pi,P2, ■ ■ ■ )Pn} be a set of n points in the plane. The /3-skeleton of P 
is a geometric graph (P, E) such that for every pi,pj G P, the edge pPpJ G P if and 
only if the /3-lune lune0{pi,pj) is empty. 

We discuss a sweepline algorithm for constructing a /3-skeleton for any /3 > 
0. The time-complexity of our algorithm is in log n) for > 1 and in 

log n), for ^ G [0, 1) under the metric Lp for 1 < p < co. In the metrics L\ 
and Loo, our algorithm takes log n) time for any /3 > 0. 
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The concept of a /3-skelcton have natural generalizations to that of A:/?-skeletons 
and additively weighted /5-skeletons. We have shown that the above sweepline 
method c&n be extended to construct these gc^ometric graphs too. The time com- 
plexity of computing the A;,b-skeleton is in 0{krv‘l‘^\og n + k?n) for /? > 1 and in 
log n 4- n^k) for (5 in the range [0, 1). The time complexity of computing the 
weighted /5-skeletoii is in log n) for /i > 1 and in log n) for /3 in the 

range [0, 1). 

We have obtained an optimal output-sensitive algorithm for computing a /5- 
skeleton for any /3 > 2 in metric Li and Lqo- The time complexity is in 0{n log n -1- 
fc), where k is the size of the /5-skeleton. 

We have presented efficient algorithms for computing the entire spectrum of fi- 
shapes in L 2 metric. This means computing for each pair of points the largest 
/5 value, Pmax, for which the /5-neighborhood is empty. Our algorithms are in 
Q(„:V 2 p^g + K) for Prnar valucs uot Icss than 1 and in 0(n^ log n + K) time 

for values greater than or equal to 0. A point that lies in the yd-neighborhood of an 
edg(! is a witness to this edge. The quantity K is a count of the number of times that 
the edg(!S of tin; initial graph on P are witnessed for /5 = 00 . We have also presented 
algorithms for computing tin; spectra of of A:/5-shapes and additively-weighted /5- 
shapes. 

We have <^xtend the definition of /5-skeleton into higher-dimensions. We have 
shown that th(^ size of th(; /5-skeleton for /5 > 2 is linear in d-dimensional space, for 
any fixed d. We have presented an efficient algorithm for computing the /5-skeletons 
for /5 > 2. 
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Chapter 1 


Introduction 


Shape is only operationally defined. Thus things do not have a shape the 
way Santa Claus has a red suit - Jan Koenderink in Solid Shape [Koe90] 


In many applications in image processing, computer vision, pattern recognition, 
computational morphology, and geometric modeling et al, an object is specified in 
terms of a set of points. The shape reconstruction and subsequent identification of 
this object from the set of input points turns out to be one of the main problems 
confronting researchers working in these areas. Similarly, in spatial analysis [GB78, 
HCF77], the locational identifiers of some real spatial phenomena are represented 
by a set of points. Given the input points, the problem here is to understand and 
characterize the real phenomena. 

A set of points, however, is an ambiguous representation of an object or a 
real phenomenon and, therefore, cannot be used directly in many applications. 
On the other hand, though the whole boundary can be used to define an object 
unambiguously, such a description would take up a huge amount of storage and a 
lot of processing time. To optimize the use of resources, we therefore need efficient 
methods to extract the object shape from the point set and fast methods to perform 
various operations on the reconstructed object. 


1 



2 


The complexity of reconstructing the shape of an object from a set of points also 
depends on the point extraction method. For example, the problem becomes harder 
if the input points are unorganized. In many applications, the point extraction 
method is not known and input points are unorganized. So an unified method for 
reconstructing the object shape is required l.o handle data from various sources. 

In some applications, the input primitives are more complex then points, through 
not unorganized. For example, a well-known problem in the field of Medical Imaging 
is to reconstruct a surface from a set of parallel cross-sectional data [BS94, GKDM96, 
AJM97]. 

It is known that, in many applications geometric properties alone cannot de- 
termine the shape of an object [TouSOb]. Even the exact solution of the geometric 
prol)lcms gives only partial, approximate, or heuristic solutions to various appli- 
cation problems. Shape reconst.ruction using non-geometric features of points is 
beyond the scope of the thesis. 

The problem of shape, recovery that we have discussed so far gives rise to the 
closely-associat(;d question of wliat if anything is meant by the shape of a set of 
points. The question is important because the tools available for the shape recovery 
problem that we have discuss(xl so far deiKuids crucially on our answer to this 
qucjstion. 

There is no single definitive answer. One reason for this is that there is no one 
definition that is best for all applications. So various definitions are extant in the 
literature. For example, Kirkpatrick and Radke [KR85] have defined the shape of 
a set of points as the decomposition of a point set into components that identify 
“essential” points and, among these, join “essential” neighbors. As is apparent, this 
definition depends on the interpretation of the term “essential” . 

Shape descriptors are divided into two categories: those that capture, so to say, 
the “external” shape of a set of points as opposed to those that capture what is 
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called the “internal” shape [KR85]. We explore these definitions in more details in 
the following sections. 


1.1 External Shapes 

The external shape of a set of points is obtained by identifying the “es- 
sential” extreme points of the set and, among these, joining “essential” 
neighbors [KR85]. 


Among the simplest structures that capture the external shape of a planar set 
of points are the minimum bounding box and the minimum spanning circle 
of this set. The minimum bounding box is the smallest rectangle which is parallel 
to the axies and encloses all the points. This can be computed in 0{n) time. The 
minimum spanning circle is the smallest circle that encloses all the points. This 
can be computed in linear time [Meg83]. Instead of a minimum spanning circle 
one could consider as an alternative the minimum spanning ellipse. This can be 
computed in O(n^) time [Pos81, Pos84]. 

In the direction of increasing descriptive complexity and refinement, next to the 
above shape descriptors is the convex hull of the given point set. The convex 
hull is the smallest convex set that encloses all the points. Computing the convex 
hull of a point set is one of the most well-researched problems in computational 
geometry [PS85]. Various properties of the convex hull are used in pattern recogni- 
tion applications (see survey [Tou80b]). Depending on the application, the convex 
hull may either give too detailed a shape description or may not give enough. 
Edelsbrunner et al [EKS83, EM94] overcame this lacuna by a natural generalization 
of convex hull to the notion of a-shapes. Here a is a real-valued parameter; so we 
get a spectrum of shapes by varying a. 

The o-shape of a point set is defined in terms of generalized discs. An a-disc, 
for an arbitrary real value a, is defined as follows: 



• if a > 0, it is a closed disc of radius l/o:; 

• if a = 0, it is a closed half-plane; and 

• if a < 0, it is the closed complement of a disc of radius —l/o:. 

For a given a, the a-hull of a point set P is the intersection of all a-discs that 
contain P. For o; = 0, the a-hull is exactly the convex hull. An a-shape of a point 
set is the boundary of the Q;-hull, with curved edges being replaced by straight edges. 
As shown in Fig. 1, a large value of a tends to produce a crude shape of the points, 
whereas a smaller value of a gives a more detailed shape. 


lA ' 

¥ 

Figure 1: a-shape for various valuesm of a. (Taken from [Ede87].) 

The notion of a-shapes has been further generalized to weighted a-shapes 
[Ede92, VBW94]. 

An a-shape is inadequate for a point set of variable density. For example, when 
we have two clusters of different densities. Melkemi [Mel97] proposed the notion of 
A-shapes to address this situation. The A-shape of a set of points P is defined by 
means of a set A of control points as follows. 

Let i?(A,pi) be the Voronoi region of pi G PU A. Let A' denote the extremities 
of the Voronoi edges that separate points in set P from those in A in the Voronoi 
diagram Vor{P U A). For q in A', let Cdisk{q) be the closed complement of the 
Delaunay disc centered at q. 
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The ^-hull of P, n{A,P) = [n,e^,CcJz5A:(g)] n B{A,P), where B{A,P) = 
\j,^sR{Ap) 

A Delaunay edge piPj G DT{PuA) is an ^-edge if there exists a Delaunay circle 
passing through pi, pj, and a point in A. The Ashape of P is the set of hedges. 
See Fig. 2. The tinne required to compute an A-shape is in 0(n log n) [Mel97]. 




Figure 2; A-shape of a set of points. (Tahen from [Mel97].) 


1.2 Internal Shapes 

The internal shape of a set of points is the shape exhibited by identify- 
ing the “essential” internal points of the set and, among these, joining 
“essential” neighbors [KR85]. 

The most elementary internal shape descriptor is obtained by joining all closest 
pairs of points. The closest pairs of points give an approximate density of the point 
set. This can be computed in d{n log n) time [PS85]. Another elementary shape 
descriptor is the graph obtained by joining each point to all its nearest neighbors. 
This graph is called the Nearest Neighbor Graph (NNG). The NNG of a point set 
can be computed in 0(n log n) time [PS85]. 

At the next level of complexity, some well-known internal shape descriptors are: 
The Sphere of Influence Graph (SIG) [Tou88], The Relative Neighborhood 
Graph (RNG) [Tou80c], The Gabriel Graph (GG) [MS84], and The Delaunay 
Triangulation (DT) [Del34]. We discuss these in detail below. 
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For each point pj, let r-j be the distance to a closest point. The circle of influence 
of Pi is the circle of radius Ti, centered at pi. The Sphere-of-Influence graph is 
obtained by joining each pair of points pi and Pj whose circles of influence intersect 
in more than one point. This can be computed in 0{n log n) time [Tou88]. 

Two points Pi and pj are rehitive neighbors if and only if the following inequality 
holds. d{pi,pj) < max{d{pi,pk),d(j}j,pk)\pk ^ P — {pi^Pj}}. Intuitively, this means 
that there is no point that is closer to both pi and p,- than they are from each 
other. The Relative Neighborhood Graph (RNG) is a geometric graph obtained by 
joining each pair of points that are relative neighbors. The RNG can be computed 
in 0{n log n) time [Sup83, HwaOO] The RNG is a subgraph of the Delaunay trian- 
gulation (DT). Given the DT of a point set, the RNG can be computed in linear 
time [Lin94, JK87]. 

Two points Pi and pj are Gabriel neighbors if and only if the circle of radius 
d{pi,pj)/2 passing through the points pi and pj is empty. The graph obtained by 
joining all Gabriel neighbors is called the Gabriel Graph (GG). It can be computed 
in 0{n log n) time [MS84]. The GG can be computed in linear time if the Delaunay 
triangulation of the point set is given [JT92]. 

The Delaunay Triangulation (DT) is a maximal planar descriptor of internal 
structure. It, has several interesting properties [PS85]. Two points pi and pj are 
Delaunay mughbors if there exists an empty circle passing through them. The 
Delaunay triangulation is a graph obtained by joining all Delaunay neighbors. The 
scope of applications of this geometric gra])h embrace fields as diverse as biol- 
ogy, crystallography, (?cology, to name a few. It can be computed in 0{n log n) 
time [PS85]. 

Just as the a-shapes give rise to a continuous spectrum of external shapes by 
varying the parameter o;, /J-skeletons, introduced by [KR85], give rise to a continuous 
spectrum of internal shapes by varying the real parameter /3. This can be used to 
extract the boundary of an objects [ABE97]. 
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More specifically, a ^-skch'ton is a geometric graph, obtained by joining pairs 
ol points whose /i-neighborlK)ods are empty. The neighborhood size depends on the 
parameter /?, and are of two types. 

The circle-based /!l-neighl)orhood of a pair of points {pi, pj) is defined as follows: 
for d > 1, it is the union of all o|)en circles of radius P\p^\/2, passing through pi 
and Pj. For /? € [0, 1], it is the interior of the intersection of two open circles of 
radius \piPj\/ (2/?), passing through pi and pj. See Fig. 3. For /? = 1 is exactly equal 
to the Gabriel graph. Given the DT of the point set, the circle based /^-skeleton can 
be computed in linear time [KR,85], for /? > 1, 

13 = 2 
(3 = 1.5 

P = 1 
(3 = 0.5 


Figure 3: Circhi-ljascd /3-neighbor hood (pi,pj) for various (i > 0. 

The lune-based /3-neighborhood of a pair {puPj) of points is defined as follows: 
for /? > 1, it is the interior of the intersection of two circles of radius ^lpjp7l/2, 
centered at the points (1 — /3/2)7>i + {0/2)pj and {p/2)pi + (1 — /3/2)pj, respec- 
tively. For P e [0, 1], the lune-based neighborhood is identical with th(j circle-based 
mughborhood. See Fig. 4. 
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P = oo 
P = 3 
P = 2 

P = l 


P = 0.75 


Figure 4; Lune-based /3-neighborhood(pi,py) for various /3 > 0. 


For ^ = 1 is exactly equal to the Gabriel graph and for /? = 2 it coincides with 
the RNG. The /5-skeleton for 1 < /5 < 2 can be computed from the DT in linear 
time [Lin94, JKY] in the metric Lp, provided p = 2 or /5 = 2. 


1.2.1 Other Shape Descriptors 

Another important shape descriptor which captures internal as well as external shape 
of a point set is the 7-skeleton [Vel92b, Vel92a, Vel94]. It is defined in terms of two 
parameters cq and Ci, where Co,ci e [-1,1] and |co| < lci|. The 7-neighborhood of 
a pair of points pi and pj is defined by two circles of radius dipi,pj)/2{l - [col), and 
d{pi,pj)l2{l — Icil) passing through pi and pj such that 

• if CqCi < 0 : the centers of the circles lie on the same side of the line passing 
through Pi,Pj, 
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• if Co Cl > 0 : the centers of the circles lie on opposite sides of the line passing 
through pi,pj, 

• if Cl € [—1,0] : the 7-neighborhood is the interior of the intersection of these 
circles, plus the interior of the edge connecting co with Ci, 

• if Cl G [0, 1] : the 7-neighborhood is the interior of the union of these circles, 
plus the interior of the edge connecting cq with Ci, 


For given points Pi,Pj and Co,Ci satisfying the above conditions, there are two 
7-neighborhoods. Two points are 7-neighbors if and only if at least one of the 
7-neighborhood is empty. The 7-skeleton is a graph obtained by joining all the 
7-neighbors. See Fig. 5 . 



Figure 5: 7(-0.15,0.3)-skeleton of a set of points. (Taken from [Vel94].) 


1.3 Organization of Thesis 


In this thesis we study some of the algorithmic aspects of lune-based /?-skeletons- 
In the next chapter, we discuss a sweepline algorithm for constructing a / 5 -skeleton 
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for any P > 0. The time-complexity of our algorithms is in log n) for /? > 1 

and in log n), for P e [0,1) under the metric Lp for 1 < p < oo. In the 

metrics Li and Loo, our algorithm takes log n) time for any P > 1. The 

concept of a /3-skeleton have natural generalizations to that of /^y^-skeletons and 
additively weighted /3-skeletons. We have shown that the above sweepline method 
can be extended to construct these geometric graphs too. In chapter 3 we present 
two output-sensitive algorithms for computing a /3-skeleton for /3 > 2, in the metries 
Lx and 33 qq. 

In the fourth chapter we present efficient algorithms for computing the entire 
spectrum of /3-shapes. Our algorithms are in log n -T LT) for /3 > 1 and in 

0{rP log n+K) time for ,5 > 0 in L 2 metric. A point that lies in the /3-neighborhood 
of an edge is a witness to this edge. The quantity K isa. count of the number of times 
that the edges of the initial graph on P are witnessed for /3 = 00. We also present 
algorithms for computing the spectra of of A:/3-shapes and additively-weighted /3- 
shapes. 

In chapter 5, we extend the definition of /3-skeleton to higher-dimensions. We 
show that the size of a ^-skeleton for ^ > 2 in higher dimensions is linear. We also 
present an efficient algorithm for computing a /3-skeleton for P > 2. 


We conclude and discuss directions for further research in the last chapter. 



Chapter 2 


Computing /5-Skeletons and Their 
Relatives 


2.1 Introduction 


In this chapter, we present a sweepline algorithm [PS85] for constructing a /5-skeleton 
for any /5 > 0. The time-complexity of our algorithm is in log n) for /5 > 1 

and in log n), for /5 G [0, 1) under the metric Xp for 1 < p < oo. In the 

metrics Li and Loo, our algorithm takes log n) time for any /5 > 0. 

The concept of a ^-skeleton have natural generalizations to that of fc^-skeletons 
and additively weighted /5-skeletons. We show that the above method can be 
extended to construct these geometric graphs too. 

The above algorithms are based on an efficient solution to the following n-circle 
problem: Given a set (7 of n circles and a set P of n points, report all the circles 
in C that are empty. 

This chapter is organized as follows. In the next section, we present an efficient 
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algorithm for the ?7.-circle proljle.m. In the third section we show how this can 
])(' used to compute a /^-skeleton for any p > 0. In the following sections, we 
define the notions of a A;/?-skeleton and additively weighted /3-skeleton. Algorithms 
for computing these are presented in the fourth and fifth sections respectively; we 
conclude in the sixth section. 


2.2 The n-Circle Problem 

2.2.1 Preliminaries 

The fi-circle problem we stated in the introduction can be solved as follows: Compute 
the Voronoi diagram of the query points. Let Vi denote the Voronoi region of a point 
p,. For each circle Cj, we find the Voronoi region Vi that contains the center of Cj. 
T1k> circle Cj is non-empty, if the distance between pi and its center is less than its 
radius, that is pj G Cj. The time complexity of this algorithm is in 0(n log n). But 
f,h(' i)rol)lem with this metliod is that checking the emptiness of limes is no longer 
as simple. 

Another method of solving the n-circle problem is by repeated point location 
<}U(!ry in an arrangement of n, circles. We could use Bentley and Maurer’s algorithm 
for tiiis purpose [BM79], which preprocesses the circles in O(n^) time for point 
location. It takes a query point, and outputs the circles that contains this point in 
0(log n -f K) time, where K is size of the output. 

For the problem at hand, however, we know all the query points off-line. Can 
we make use of this fact to solve our problem more efficiently? The following 
observations show that the answer is yes. 

Observation 2.1 An arrangement of n circles divides the plane into a set R of at 
most (){ 7 p) regions. 
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Proof: Let n - 1 circles divide the plane into T{n - 1) regions. The circle 

intersects at most n - 1 circles. These intersection points divide the boundary of 
the circle into at most 2(n - 1) disjoint arcs, since there are at most 2(n - 1) 
intersection points and two consecutive intersection points define an arc. Each arc 
of the circle divides a region into two. So the extra regions created by the 
circle is at most 2(n — 1). Therefore T(n) = r(n - 1) + 0(2n - 2). This implies 
that T(n) = O(n^), since T(l)= 1. | 


Observation 2.2 If a circle is not empty, then it is sufficient to find one witness 
(that is, a point lying inside the circle) for it. 


Observation 2.3 If a point is in the region r £ R, then all the circles that contain 
this region are non-empty and can be pruned. 


2.2.2 An Algorithm 

The main idea of our algorithm is that instead of pre-computing the arrangement and 
then doing repeated point location, we can achieve the same effect more efficiently 
by sweeping the circles and the points simultaneously. For details on the sweepline 
paradigm we refer the reader to [PS85]. 

We shall use the following notations and definitions throughout this chapter. 
For each position of the sweepline, we define a partial order relation <„ in the set 
of regions R as follows. For two regions n and rj in R, ri <a rj if the intersection of 
the region rj with vertical line x = a lies below the intersection of the region rj with 
that line. This ordering information is maintained in the sweepline status C. An 
event point is either an intersection point of two circles, or the leftmost/rightmost 
point of a circle, or an input point. The events are maintained in a priority queue Q, 
which is initialized to contain the input points and the left-end points of the circles. 
A circle in C is said to be active if it is either empty or a witness is yet to be found. 
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We denote the set of active circles by Ac. Initially, Ac = C. For each region r e C, 
let Cr denote the set of active circles that contain this region. 

To describe the algorithm it is enough to specify the actions at each type of 
event point as the sweep proceeds. 

Let the event e correspond to the left end point of a circle c. We insert the 
right-end point of the circle in Q and locate the region r e C that contains this 
event point. The boundary of c splits r into three regions ri, r 2 , and r^, as shown 
in Fig. 6. 

Comment : It may be noted in passing that ri and are the same region r as can 
also be seen from the Fig. 6; however, we distinguish between them because in the 
sweepline status a region is identified and maintained by the two circular arcs that 
bound it from above and below. Thus r appears in two incarnations in the sweepline 
status in this case. i 

Therefore, we replace the region r by the regions ri, r 2 , and r^ in £. Simultane- 
ously, we create the active circle lists Cm Cm and Cr^\ the first and third are really 
copies of Cr, while the second is obtained from Cr by adding to it the circle just 
encountered. If c intersects the circular arc that bounds ri from above we insert this 
intersection point in Q. Similarly, for the circular arc that bounds tz from below. 
These intersection points are qi and in Fig. 6. 

Let the event e corresponds to the right end point of a circle c. If this point is 
on the boundary of the regions ri, r 2 , and rz, then the region r 2 is deleted from £ 
and the regions ri and rz are coalesced into a new region riz (see Fig. 7). The active 
circle list C^ of this new region is the same as that of ri or rz. We insert into 
the event queue the intersection point, if any, of the bounding arcs of the region riz 
(point q in Fig. 7). This intersection point can very well be the right-end point of a 
circle as shown in Fig. 8. 

Let the event e correspond to an intersection point of two circles Ci and C 2 . If e 
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is on the boundary of the regions ri, r^, and r^, then these regions transform into 
new regions r[, r'^ and respectively. The region r'l has a new lower boundary 
vis-a-vis and the same active circle list; has new upper and lower boundaries 
as compared to and an active circle list obtained from Crj by deleting/adding 
the intersecting circles C\ and C 2 ; rg has a new upper boundary compared with 
and the same active circle list. If the bounding arcs of r[ intersect, we insert their 
intersection point in Q. We do the same for rg. In Fig. 9 these intersection points 
are qi and 92 - Notice that the latter is the right end point of one of the intersecting 
circles. 



Figure 9: Update of C at an intersection point of two circles. 


Let the event e correspond to an input point, we locate the region t ^ C, that 
contains this point and remove all the circles in Cr (see Fig. 10). The removal of 
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each such circle involves deleting it from the set of active circles Ac, deleting the 
intersection points of circle c with other circles, which are present in Q, and updating 
C. 



Figure 10; Update of C at an input point. 


2.2.3 Analysis of the Algorithm 

We claim the following result. 

Lemma 2.1 The n-circle problem can he solved in 0{n^ log n) time, using O(n^) 
space. 

Proof: At any abscissa, the sweepline has at most 2n intersections with the 

circles; hence at most 2n + 1 segments. Since each segment corresponds to a region, 
the sweepline status can have at most 0(n) regions. 

Any data structure used for implementing the sweepline status £ should support 
efficient implementations of the following operations: addition and deletion of a 
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legion, point location, and merging of two adjacent regions. A height-balanced tree 
meets these needs eminently. The regions are stored at the leaves and the internal 
nodes correspond to the circular arcs' that make-up the region boundaries. The 
active ciicle list Cr of each region r in the sweepline status is also maintained as a 
height-balanced tree. 

We maintain a Boolean array Pruned to keep track of the status (empty /non- 
empty) of all the circles. Initially, Prune[i\ is false for each i in 1 < f < n. During 
the sweep, if circle Cj is found to be non-empty we set Pr'une[i] to true. 

The event queue, as well as the array Prune[n], can be initialized in 0(n) time. 

We analyze separately the complexity of processing each of the four types of 
events in Q. 

Let e be an event that is the left end point of a circle c. The cost of locating 
this point in C is in 0(log n). If r be the region that contains this event point, 
we replace it by three new regions ri, r^, and vz in £ and if the bounding arcs of 
the region ri intersect, insert their intersection point in Q. We do the same for the 
region r^. The complexity of this is in 0(log n) again. Creating the active circle lists 
corresponding to these regions is in 0(n), since we have to make three copies of Cr 
and insert the newly encountered circle into one of these. The cost of processing the 
left-end point of a circle is in 0{n) time. The total cost of processing the left-end 
points of all the circles is in O(n^) time. 

A similar analysis shows that the cost of processing the right-end point of a 
circle is in 0(log n). The total cost of processing the right-end points of all the 
circles is in 0(n log n). 

Let the event e be an intersection point. The cost of locating this point in £ and 
hence the three consecutive regions, ri, r 2 , and rs, that have e is on their boundary is 
in 0(log n). The cost of replacing these regions in £ by the corresponding regions r(, 
Tj and is also in 0(log n). The cost of creating the active circle lists corresponding 
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to these regions is in 0(log n), since we only need to make two deletions/insertions 
into the list corresponding to rz to get the active circle list of r'^. The active circle 
lists of r[ and r'^ are same as that of ri and respectively. We insert in Q, the 
intersection points of the arcs that bound the region r[ as discussed earlier. We do 
the same for r^. The complexity of this is in 0(log n) again. So the cost of processing 
one intersection point is in 0(log n) time. The cost, therefore, of processing O(n^) 
intersection points is in 0(n^ log n). 

Let the event e be an input point. The complexity of locating this point in the 
sweepline status is in 0(log n). If r be the region that includes this point, we delete 
all the circles from its active list Cr and set the appropriate array entries in Prune\} 
to true. 

If the point e is on the boundary of the circle c and one of the bounding arc of 
the region r is this circle c. Let the regions ri and r 2 are respectively above and 
below the point e in the sweepline status C (see Fig. 11). The circle lists of these 
regions are differ by the circle c. That is, — {c} = Crj or Crj — {c} = ■ The 

point e is not a witness for c, but it is a witness of all other circles of the list 
and Crj. Therefore, we prune all the circles in — {c}. 



Figure 11: Update of C at an input point, if it is on the boundary of a region. 
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The cost of deleting a circle Ci is in 0{'n log 7 i). We argue as follows. As 
observed earlier, I his involves removing from the event, queue its intersections with 
othei ciicles that, aie present in Q as also its right end-point; it also involves removing 
this circle from tin; active list of all the regions in £ in which it occurs. The first can 
be done in constant time by the following artifice. We merely set Prune[i] to true. 
While processing an intersection point g or a right end point we have to take the 
additional care of checking suitable arraj' entries in the array Pruned to decide if 
q needs to be processed. The cost of the second is in 0(ri log n). The adjustments 
required to £ is in 0(log n) since the removal of its upper /lower arc causes the 
regions above and below this arc to be merged into one. In the course of the sweep 
at most n circh's are pruned, and thus the amortized cost is in 0{n^ log n). 

Hence the n-circle problem can be solved in 0(n^ log n) time using O(n^) space. 

I 


The comjdexity of the above algorithm is worse than the brute-force algorithm 
which tests a i)oiiit. for inclusion in each of the n circles! For the complexity of this 
is in 0(n), and tluhefore in 0{iP) for n points. However, if we look at the above 
algorithm very canifully we notice that the problem lies in considering all the circles 
in one go. Tlu' processing time for the n circles is high and dominates the algorithm. 
The following lemma shows how to balance the processing time for circles and points 
to imj)rov(^ the time’ complexity. 


Lemma 2.2 Given a .set C of n circles and a set P of n points, all the empty circles 
can be comptUed in OipifG log 7^) time using 0(n) space. 


Proof: We divide the problem into n/m subproblems, each subproblem consisting 

of rn circles and al! n input' points. For each of these subproblems the above 
algorithm takes (){nif log m + n log m) time. Therefore, the total time required 
to solve all the n/m ])roblems is 0(nm log m -f n^ jm. log m) time. The minimum 
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time of log n) is achieved for m = -y/n. The space required is in 0{n), since 

each subproblem consists of circles. Hence the claim of the lemma. i 


2.3 /3-Skeleton 


In this section, we generalize the above algorithm to compute the /5-skeleton for 

/?>!. 

It is well-known that the Delaunay triangulation (DT) [PS85], is a super-graph 
of a /5-skeleton for /5 > 1. Therefore, to compute the /5-skeleton, for some /3>1, we 
construct DT(P) and prune the edges whose yd-lunes are not empty. The Delaunay 
triangulation can be constructed by any standard algorithm in 0(n log n) time. 
What remains is to find a way to prune the non-empty lunes. 

It is easy to extend the solution to the n-circle problem, if circles are replaced by 
lunes. The same observations have been made by Su and Chang in [SC91b]. In fact, 
we can consider any set of convex objects with constant-size description, any two of 
which have a bounded number of intersections. With lunes, the event structure is 
slightly more complicated. Unlike in the case of circles, each lune can have up to 
four events associated with it. Two of these correspond to positions where a lune 
joins the sweepline and leaves it. The other two correspond to the intersections of 
the circles that define a lune (ii and in Fig. 12). 

Due to this event structure the region in the sweepline status that reflects its 
membership in £ can be defined by three different pairs of arcs. For example in 
Fig. 12 these are the pair {a, 6}, {a,c} and {c,d} consecutively. The rest of the 
details carry over and we have the following claim: 

Theorem 2.3 For /5 > 1, the (3 -skeleton of P can be constructed in 0{n^^^ log n) 
time using 0{n) space. 
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Figure 12: Initial events of a lune. 


Similarly, the /?-skeleton of P, for /? € [0, 1), can be constructed by starting with 
the complete graph on P as the super graph and pruning the non-empty lunes. The 
following theorem can be proved along exactly the same lines 


Theorem 2.4 The 13-skeleton of a set P of n points in the plane, for P e [0,1), 
can be constructed in 0{n^^‘^ log n) time using 0{n^) space. 


Proof: We divide Oin^) lunes corresponding to the edges in the complete graph 

into n^/m subproblems. Each subproblem consists of m lunes and all n input 
points. For each of these subproblems the above algorithm takes 0{nrp log m + 
n log m) time. Therefore, the total time required to solve all the rP/m problems is 
0{n^m log m + nVm log m) time. The minimum time of log n) is achieved 

for m = y/n. The space required is in 0{n^), since we are start with complete graph. 
Hence the claim of the lemma. I 


2.4 A;/3-Skeleton 

We can generalize the notion of a /^-skeleton by relaxing the emptiness criterion and 
allow the interior of a lune to contain up to A: - 1 points; we get what is called a 
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A', 5-skeleton. The geometric graphs that we get for different values of 1 upwards, 
are all subgraphs of what is called the A-Delaunay graph. 

The A-Delaunay graph is a generalization of the Delaunay triangulation on n 
points. The triangulation in this case can be characterized this way: a triangle 
belongs to the A-Delaunay graph if its circumcircle contains at most A - 1 points. 
^\e can obtain it from the Ath-order Voronoi diagram as follows: an edge PiPj belongs 
to the A-Delaunay graph if a segment of the bisector of pi and pj is in the Ath-order 
Voronoi diagram [SC90]. 

Algorithms for computing A/?-skeletons for special values of p are extant in the 
literature. Su and Chang [SC90] has shown how to compute the A-Gabriel graph 
(,5 = 1) in C)(A^ n log n) time. For fixed values of A, Chang et al. [CTL92b] presented 
an O(n^) algorithm to compute the A-Relative Neighborhood Graph (/? = 2). Later, 
Su et al. improved this to log n) [SC91b]. An interesting application of the 

A-RNG has been made to the Euclidean bottleneck matching problem, by Chang et 
al [CTL92b]. In this problem a matching has to be found among n points in the 
Euclidean plane that minimizes the longest edge. 

To compute a A/3-skeleton for /3 > 1, we start with the A-Delaunay graph of P 
and prune those edges whose lunes contain more than A — 1 points. It turns out 
that the above n-circle algorithm can be adapted to this case also. The additional 
work involved is to maintain a counter for the number of points inside each lune. 
We increment this count for the appropriate set of lunes whenever an input point 
is encountered during the sweep. Whenever the counter value of any lune becomes 
A, we prune that lune from further processing. The extra time required to maintain 
all the counters is in 0{k^n), since there are 0(A(n - A)) counters (the A-Delaunay 
graph has 0{k{n - k)) edges [SC90]) and each counter is incremented at most A 
times. We have thus proved the following theorem. 

Theorem 2.5 The kp-skeleton of a set P of n points in the plane for p>l can be 
constructed in 0{kn^l^ log n + A^n) time. 
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It is quite easy to compute a A:/?-skeleton, for values of /3 in the range [0,1). We 
need to start with the complete graph on P, and prune edges which are not in the 
/c/3-skeleton. The additional cost of this is reflected in the following result. 


Theorem 2.6 The kp-skeleton of a set P of n points in the plane, for /? 6 [0, 1) 
can be constructed in log n + n'^k) time. 


2.5 Weighted Neighborhood Graphs 


In this section, we discuss various additively weighted neighborhood graphs. We 
assign a positive weight Wi to each point pi e P. The weight Wi expresses the power 
of the point pi in the sense that each weighted point pi can be treated as a circle 
Bi, centered at pi, with radius equal to the weight Wi. Thus the weighted distance 
dw{Pi,^), between a point Pi G P and an arbitrary point x, is d(pi,x) — Wi if x is 
outside Si, and zero otherwise. The weighted distance dyj{pi,pj), between any two 
points Pi,Pj € P is the distance between Sj and Bj. Let D denote the set of circles 
Bi. We discuss the details of a few individual special cases 

2.5.1 Weighted Gabriel Graph (WGG) 

The weighted Gabriel neighborhood, NoipuPj), of two points pi,pj e P that 
have disjoint circles S, and Bj is an open circle of radius dy,{pi,pj)/2, whose center 
is on the line segment Wj and is tangent to the circles Bi and Bj (see Fig. 13). If 
Bi and Bj intersect, then NaiPuPj) is empty. 

The weighted Gabriel graph of a set of points P is a geometric graph (P, E) 
such that for every puPj 6 P, the edge piPj € P if and only if No{pi,Pj) is not 
intersected by any circle Bk G D — {Bi,Bj}. 
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The weighted Delaunay graph is a dual of the weighted Voronoi diagram [AE84] 
and can be computed in 0(n log n) time [Mir93]. An edge belongs to this 
graph, if there is a common boundary between the Voronoi regions of pi and pj. 
The following result and its proof are from Mirzaian [Mir93]. 


Lemma 2.7 The weighted Gabriel graph is a subgraph of the weighted Delaunay 
graph if the set of circles D is pairwise disjoint. 


Proof: Let the edge pip] G WGG. The Gabriel neighborhood of Pi and pj, 

NciPiiPj), is not intersected by any circle in D. Let the point Cp^p^ be the center of 
the Gabriel Neighborhood Na{pi,pj). The circles Bi and Bj are at an equal distance 
from Cpipj. Moreover, no circle Bk £ D — {Bi,Bj} is closer to Cp^pj than the circles 
Bi and Bj. Prom the definition of a Voronoi diagram, the point is on the 
bisector of the sites Bi and Bj. Thus by definition p^ G WDG [Mir93]. | 


2.5.2 Weighted Relative Neighborhood Graph (WRNG) 

The weighted relative neighborhood, of two points G P that 

have disjoint circles Bi and Bj is the interior of the intersection of two circles, of 
radii dyjipuPj) + Wi and dyjipuPj) + Wj, centered respectively at the points pi and 
Pj (see Fig. 14). If Bi and Bj intersect then NRNciPuPj) is empty. 
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The weighted relative neighborhood graph of P is defined as follows. For every 
Pi,Pj € P, the edge pipj is in the graph if and only if the following inequality holds. 
dwiPi,Pj) < max{dyj{pi,pk),d^{pj,pk)}, for all p* € P - {pi,Pj} (see Fig. 14). Or 
equivalently, if the lune is not intersected by any circle in D — {Bi, Bj}. 

Mirzaian [Mir93] made a very interesting application of the weighted RNG to 
the problem of Minimum Weight Euclidean Matching. In this problem we are given 
2n points in the Euclidean plane and are required to match them into n pairs such 
that the sum of the n distances between the matched pairs is minimum. 


2.5.3 Weighted /^-Skeleton 

For all pi,pj e P such that Bi and Bj are disjoint, the weighted lune-based neigh- 
borhood, lune^{pi,pj), is a region lying between the circles Bi and Bj and is defined 
as follows: for /? > 1, it is the interior of the intersection of two circles of radii 
lddvj{pi,pj)/2 -b Wi{l3 - 1) and l3d^{puPj)/2 -b Wj{/3 - 1), whose centers are on the 
line joining Pi and Pj, and which touch externally the circles Bj and Bi respectively 
(see Fig. 15). For /? € [0,1], it is the interior of the intersection of two circles of 
radius dyj(pi,pj)/2p, passing through p[ and p'- which are a pair of closest points 
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on the circles Bi and Bj (see Fig. 15). If the circles Bi and Bj intersect, then the 
htne0{pi,pj) is empty for all P >0. 

^ /? = 3 

P = 2 

/3 = 1 

Pj 

/3 = 0.75 


Figure 15: Neighborhood of weighted /3-skeleton. 

The weighted /3-skeleton of P in this case is a geometric graph (P, E) such that 
for every PitPj € P, the edge pip] G P if and only if lune0{pi,pj) is not intersected 
by any circle Ck € D - {Bi, Bj}, where D is the set of circles Bi. 

When /3 = 1, luneffipuPj) is exactly the Gabriel neighborhood NaipuPj), and 
when ^ = 2, we get the weighted relative neighborhood, NRNoiPi-iPj)- 

Lemma 2.8 For a pairwise disjoint set of circles D and /3 > 1, the weighted 
skeleton is a subset of the WDG. 



Proof; This lemma follows from the fact that for 13 = 1, the weighted /3-skeleton 
is exactly the weighted Gabriel graph, which is a subgraph of WDG from Lemma 2.7. 
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2.5.4 Algorithm 

In this section, we discuss an algorithm for computing a weighted ,5-skeleton for 
/? > 1. We assum® that the circles in D are pairwise disjoint. 

As we mentioned in the previous section, the WDG is a super-grajih of a W/3- 
skeleton for /5 > 1. The \\ DG can be constructed in 0{n log n) time [Mir93]. 

So we just have to prune the lunes that intersect any circle Bi G D. For the 
sake of simplicity, we consider circles instead of lunes. As we observed earlier, the 
generalization to lunes is pretty straightforward. The problem can be stated with a 
bit of color as follows; 


Given a set of 0(n) blue circles C and a set of n red circles D, find all 
the blue circles that are not intersected by any red circle. 


We follow closely the basic circle-point algorithm to design an algorithm for 
this case. Most of the details carry over, except when in the course of the sweep 
we meet a red circle we locate its left-end point in a region of the sweepline status 
and check for intersection with the circles which correspond to the upper and lower 
arcs that bound this region. If a bounding arc is intersected, the corresponding 
circle is removed from the active list of blue circles, and the arc is deleted from 
the sweepline status. This latter action causes two adjacent region to be merged 
into a new one and we repeat the above sequence of steps for this region. The 
process terminates when we get a region whose bounding arcs have no intersection 
with the red circle. Now we insert the red circle into the sweepline status because 
there can be yet more intersections with blue circles that have not been swept yet. 
However, such an intersection can occur only after a blue-blue intersection has been 
encountered during the sweep (see Fig. 16). So, whenever a red-blue intersection 
is encountered, prune the blue-circle corresponding to this intersection. A routine 
complexity analysis now gives the following result. 
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Figure 16: Intersection between the red circle Bk and the blue circle Cj is detected after blue-blue 
intersection tij. 

Theorem 2.9 The weighted /3-skeleton of a set of n points in the plane for /5 > 1 
can be constructed in log n) time. 

Theorem 2.10 The weighted ^-skeleton of a set of n points in the plane for (3. € 
[0, 1) can be constructed in log n) time. 

2.6 Conclusion 

In this chapter, we have presented an efficient sweepline algorithm for the n-circle 
problem, and used this to compute /^-skeleton, A:/?-skeleton, and weighted jd-skeleton 
under the metric Lp for 1 < p < oo. In the metrics Li and L^o, the size of 
these geometric graphs are in 0{n^) for any p. We can compute these graphs in 
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0(n^/^ log n) time by starting with the complete graph and pruning the non-empty 
lunes, following the technique of this chapter. 

The algorithms presented in this chapter are not optimal. It would be intersect- 
ing to find out more efficient algorithms for computing these geometric graphs. 



Chapter 3 


An Output- Sensitive Algorithm 
for Computing the /3-Skeleton 

3.1 Introduction 


In the previous chapter we have presented an efficient algorithm for computing the 
/3-skeleton in any metric Lp, for 1 < p < oo. In the metrics Lx and Loo, a /?- 
neighborhood is a simple rectangle compared to the complex lune in metric Lp for 
1 < p < oc. Now the question is can we make use of this property to design more 
efficient algorithms? The answer to this question is yes. We present two output 
sensitive algorithms for computing a ^^-skeleton in the metrics Li and Loo for any 
/3 > 2. The algorithms are in O(n^/^log n + k) and 0{n log n + k) respectively, 
where k is size of the /3-skeleton. For the rest of this discussion we shall confine 
ourselves to the Loo metric because rotating by 45° and scaling by a factor of s/2, 
the Li metric is transformed into Loo- 
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3.2 Size of a /5-skeleton 


In the metric Lp for 1 < p < oo, a ;5-skeleton is a subgraph of the Delaunay 
triangulation [TouSOc] for p > 1. This implies that the size of the /5-skeleton is in 
0{n), for p > 1. But in the Loo metric , the size of a /3-skeleton is in Q{tP). We 
show this by constructing an example. 

Let us place n points on each vertical edge of a square as shown in Fig. 17. The /3- 
neighborhood of each edge pq, where p E {pi, , Pi^, ... ,pij &nd q e {pj, , Pj, , ■ . . , Pj„ } , 
is empty for any /3 > 1. Therefore, the size of a jd-skeleton for any /3 > 1 is in 0{rP). 



3.3 Algorithm 


The /3-skeleton can be computed in 0{rP) time by starting with the complete graph 
as a super graph and pruning the non-empty lunes by a brute force method. This can 
be improved to log n) time by the sweepline method described in previous 

chapter. 

In the metrics Li and Loo, the /3-lune is a rectangle (see Fig. 18). Therefore, the 
emptiness query of a lune can be decided efficiently in 0(log tl) time using range 
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tree data structure [PS85]. Hence a ^-skeleton can be computed in 0{v? log n) 
time. By taking a closer look into the geometry of the problem, we can improve on 
this complexity further. 


/? = 2 




/? = 3 



Figure 18: The /3-neighborhood in metric Loo for ^ > 2. 


The main idea is to find a smaller super graph than the complete graph on P. 
We adopt Yao’s [Yao82] region approach for this. 
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Consider a point pi € P. We divide the space around pj into eight regions [Yao82], 
by drawing four lines through pi that make an angle of 0°, 45°, 90°, and 135° respec- 
tively with the x-axis. These regions are open on one side and closed on the other. 
For example, Ri is closed on the x-axis and open on the diagonal, and R 2 is closed 
on the diagonal and open on the y-axis. The other regions are similarly defined. 
These regions are called narrow regions^ because for any two points Pj,Pk in a region 
dist{pj,pk) < max{dist{pi,pj),dist{pi,pk)}. Let the region, for 1 < m < 8, be 
labeled Rm{Pi) as shown in the Fig. 19. 



Figure 19: Narrow regions in metric Loo- 

We connect each point pi to all its nearest neighbors in each region RmiPi), for 
1 < m < 8. The graph so obtained is called the Geographical Neighborhood Graph 
(GNG). The following lemma is easily verified. 

Lemma 3.1 A -skeleton for j5 >2 is a subgraph of the GNG. 

Proof: Prom definition we know that any /^-skeleton for /5 > 2 is a subgraph of 

the /5-skeleton for fi = 2. So it is sufficient to show that the ^-skeleton for ^ = 2 is 
a subgraph of the GNG. 

Let assume that the edge pipj € /5-skeleton for /5 = 2. So, lune(pi,pj) is empty. 
This y3-lune is defined by two circles of radius d{puPj), centered at pi and pj. 
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Without loss of generality we assume that pj e Ri{pi). Let p^ ^ GNG. This 
implies that there exist a point p^ € Riipi) such that d{pi,pk) < d{pi,pj). That is, 
Pfc is inside the circle of radius centered at pi. 

The definition of the narrow region d{pj,pi.) < max{d{pi,pj), d{pi,pk)}, and 

d{pi,pk) < dipi,pj) implies d(p,-,p,.) < d{pi,pj). That is, pk is inside the circle of 
radius d{pi,pj), centered at pj. 

That is, pk e lunep= 2 iPi,Pj), which is a contribution. Therefore, p^ e GNG. 

Hence the claim of the lemma. | 

Now the question is, for a given point how efficiently can we find its nearest 
neighbors in each region? We consider the region Ri(pi) for a point Pi € P. The 
remaining regions P 2 , P 3 , • • • , Ps can be handled in the same way. 

We sort all the points on their x-coordinate as primary key and y-coordinate 
as secondary key. The sorted list of points is kept in an array. For each point pi, 
we find a nearest neighbor pj (in Ri{pi)) if it exists in 0(n log n) time using the 
divide and conquer approach of Guibas and Stolfi [GS83]. The rest of the nearest 
neighbors can be found as follows. 

We compute the intersection points S and N of the boundary of the region 
RiiPi) with the vertical line passing through pj (see Fig. 20). We locate the points 
S and N in the sorted array of points. All the points between N and S are the 
nearest neighbors of Pi in the region Pi(pi), since points are sorted by y-coordinate 
as the secondary key. We denote these by Pj^jPj^, ■■■,Pjk- 


3,3.1 Pruning edges 

Our algorithm for pruning the edges whose lunes are non-empty is based on the 
following key observation. 
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Figure 20: All nearest neighbors of pi in iJi. 

Observation 3.1 Let Pji be a nearest neighbor ofpi in the region R^. Then for any 
P > I, the lunep{pi,pj') is a bounded rectangle, contained within an infinite strip, 
bounded by vertical lines passing through the points pi and py . 


Let Oi and bi be points (of P) inside the strip such that lines through these, or- 
thogonal to the bounding edges of the strip, define the largest rectangle < ri, r 2 , rs, r 4 > 
contained in the strip whose interior is empty as shown in the Fig. 21. Note that 
it’s possible that one or both of these points do not exist. 

Clearly, all the lunes contained in this rectangle are empty. For /? > 2, an easy 
decision algorithm is based on the following observation. 


Observation 3.2 The regions lunep{pi,pjf), lunei 3 {pi,pjfi^, . . ., luneg{pi,pjfi) form 
a sequence of nested neighborhoods, with lune^ipuPji^i) C lunep{pi,pjf), where 1 < 
l<k-l. 

Here is the algorithm. Starting with the smallest one, luney{jpi,Pjfi), we look for 
the first non-empty lune. If none is found, all the edges between pi and its nearest 
neighbors are in the /3-skeleton. Else if lunep(pi,Pj„) is the first non-empty lune 
then the edges joining pi to Pj„ . . -Pji are pruned. 
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Obviously enough, the time required for this procedure is proportional to the 
number of non-empty lunes. 

So all that remains to be done is to show how to compute the points Uj and hi. 
Let us do this for the point bi] the other point Uj can be computed similarly. Since 
we would like to do this for every point in P, the following abstract problem can be 
formulated. 


Given a set / of n horizontal intervals and a set P of n points, for each 
interval k € /, find the closest point hi e P, if any, that lies below it 
and within the infinite strip determined by the vertical lines that pass 
through the end points of k 


Caveat: It just happens that for our problem points and intervals are not 
distinct; each point is the left end-point of an interval. 

We outline two different solutions to this problem: the first based on the familiar 
sweepline technique, and the second based on the range-tree data structure. 


3.3.2 Sweepline Method 

We sweep the plane from top to bottom. At any instant, the sweepline status £, 
contains a relative ordering of the end-points of the intervals that have been swept 
by the sweepline and whose closest point is yet to be found. Two adjacent end- 
points in C represent an elementary interval. With each elementary interval El, we 
associate a list of intervals Cei that contain this elementary interval. 

The event queue consists of the points of P in y-sorted order. 

Each stop in the event queue triggers off two events. The first corresponds to 
the point p. We locate the elementary interval El in £ that contains p. The point 
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p is closest to all the intervals in Cei- We delete these intervals from the sweepline 
status and adjust the lists of all the elementary intervals. 

The second event corresponds to the interval I whose left end-point is p. Let g 
be its right end-point. The point p divides El into two intervals Eli and El 2 - We 
delete El and insert the elementary intervals Eli and EI 2 into £. We also create 
the lists Ceii and Cei 2 - We then do the same for q. For each elementary interval 
El' that is completely inside I, we update Cep by inserting I into it. 

We claim the following result. 


Lemma 3.2 The running time of the above algorithm is in 0{n^ log n); the space 
complexity in 0(n^). 


Proof: At any time, the sweepline status £ contains at most 2n end-points, 

corresponding to n intervals. We implement £ by a height-balanced tree that stores 
the end-points at the internal nodes and the elementary intervals at the leaves of 
the tree. 

For each elementary interval El G £ we maintain in a height-balanced tree the 
set of intervals Cei that cover El. 

The event-queue is implemented by an array. It is initialized at the beginning 
of the sweep in 0{n log n) time. When the sweepline encounters an event, we first 
process the point p corresponding to this event. We locate the elementary interval 
El €C that contains p. The time for this is in a(log n). Since the point p is closest 
to all the intervals in Cei, these are deleted. If I' is one such interval, we delete 
its left-end point from £ in 0(log n) time. This causes two adjacent elementary 
intervals to be merged. The time for this is in 0(log n) . The right-end point of 
I' is similarly deleted. We also delete V from the interval list of each El, which is 
in between the left and right end points of I' in C. This can be done in 0(n log n) 
time, since there can be at most 0{n) elementary interval in sweepline. The cost of 
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deleting an interval is thus in 0[n log n) time. If we charge this cost to the interval 
itself, the total cost of deleting all the intervals is in 0{n^ log n). 

We next piocess the associated interval 1. In O(log n) time we locate the 
elementary interval El e £. that contains its left-end point a and insert it in C 
in 0(log n) time. This means splitting El into two intervals Eh and Eh, deleting 
El and inserting Eh and Eh into £. All this can be done in 0(log n) time. 

The interval list Csh is the same as Cei, while Cei^ = Cei U {1}. The cost of 
creating Csh is in 0{n). 

We repeat the above sequence of actions for the right end-point of 1. 

Further, for each elementary interval EF that is contained in I, we update Cee 
by inserting I into it. This can be done by traversing the sweepline C from the 
elementary interval, containing left-end point to the elementary interval containing 
right-end point. The cost of all these operations is in 0{n log n) time. 

Thus the the total cost of processing n intervals is in 0(n^ log n) time. 

Hence this problem can be solved in 0(n^ log n) time. | 

The complexity of the above algorithm is prima facie worse than that of the 
algorithm which, for each given interval, tests by brute-force which of the n points 
is closest to this interval. The complexity of the latter is in O(n^). However, if 
we look at the above algorithm very carefully we notice that the problem lies in 
considering all the intervals in one go. The cost of maintaining the sweepline status 
of n intervals is high and dominates the algorithm. The following lemma shows how 
the complexity can be improved by balancing the processing time for intervals and 
points. 

Lemma 3.3 Given a set I of n horizontal intervals and a set P of n points, for 
each interval li € I , the closest point h G P, which is within the vertical infinite 
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strip bounded by the lines passing through the end points of k and below k can be 
found in log n) time using 0(n) space. 


Proof: We divide the problem into njm subproblems, each subproblem consisting 

of m intervals and all n input points. For each of these subproblems the above 
algorithm takes 0(m2 log m + n log m) time. Therefore, the total time required 
to solve all the n/m subproblems is in 0(nm log m + n?' /m log m). The minimum 
time of log n) is achieved for m = ^/n. The space required is 0{n), since we 

have to consider y/n intervals at a time. Hence the claim of the lemma. | 


3.3.3 Range-tree Method 

The complexity can be further improved by using the range-tree data structure [PS85]. 
For the sake of completeness we briefly describe the method. We preprocess the set 
of points such that for a given a query interval k, the corresponding point bi, if it 
exists, can be located quickly. The range-tree is a rooted binary tree which allows 
us to search along two orthogonal directions, say x and y, simultaneously. It has 
two levels, the primary and the secondary. The primary level is useful for searching 
along the x-direction, while the secondary level for searching in the y-direction. 

Each non-root node v of the range-tree T corresponds to a set points, lying 
between the abscissae L{v) and R{v). These points are stored in an array, sorted 
by their y-coordinates, in the secondary structure Y{v). The range R{v) — L{v) 
is divided into two sub-ranges such that each sub-range contains at least half the 
number of points of v. The left and right ranges correspond to the nodes LSON{v) 
and RSON{v) respectively. The root of T corresponds to whole range, and the 
secondary level at the root node is a threaded binary tree which contains all the 
points, sorted by their y-coordinates. 
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Foi each non-leaf node v, and for each point p of K(t’), there is a pointer 
LBRIDGE{p) (respectively RBRIDGE{p)) to a point p' oiY{LSON{v)) (respec- 
tively f of y {RSON (v))), which is p itself or or a point that lies immediately below 
P- 


The range-tree can be constructed in 0{n log n) time in a bottom-up manner. 
The space required is in 0{n log n) since there are 0(log n) levels and each level 
requires 0{n) space. 

The application to our problem is quite simple. We narrow down the x-range in 
which a query interval k lies. Each time we do this, we update the point p that is 
closest to the supporting line of k in the current range to the point that is closest in 
the new range. For this latter work we use the secondary structure stored at each 
node. 

For a given query interval k, the details of the search mechanism is as fol- 
lows. Find the point p which lies immediately below the line passing through Zj by 
searching in Y{root) in O(log n) time. For each node v, if k is completely within 
the range of LSON{v) (respectively RSON{v)), we search in the sub- tree rooted 
at LSON{v) (respectively RSON{v)). Otherwise, search in both the off-springs. 
We consider searching in the left off-spring, as the right off-spring can be handled 
similarly. While proceeding to next level we can get the point which is immediately 
below the line passing through k in Y {LSON{v)), by using the pointer LBRIDGE. 
If the range of any node is completely within the given query line segment I, the 
points which are immediately below I are candidates. Of these, the required point is 
the one that is closest. See Fig. 22 and Fig. 23. The query time is thus in (9(log n) 
time. 

We summarise our results in the following lemma. 


Lemma 3.4 The algorithm for computing a fi-skeleton for P >2 in region Ri takes 
0{n log n) time. 



Figure 22: partition into standard intervals. 


The above algorithm can be applied to the remaining regions by using an 
appropriate linear transformation on the input points to get the complete /3-skeleton 
for the given /? > 2. Reporting the edges of the ^-skeleton takes time proportional 
to the size of the output. Hence follows the theorem. 

Theorem 3.5 For 0 > 2, a P-skeleton can be computed in 0{n log n + k) time, 
where k is size of the output. 


3.4 Conclusion 

In this chapter we have described an output sensitive algorithm for computing a 
/3-skeleton in the Loo metric for 0 > 2. The complexity of the algorithm is in 
0(n log n -hk). It would be interesting to find an output-sensitive algorithm for 0 
in the range [1, 2). 








Chapter 4 


Computing the /^-Spectrum and 
Their Relatives 


4.1 Introduction 


In this chapter we propose efficient algorithms for computing the entire ^-spectrum 
of a finite planar point set in the L 2 metric. This means computing for each pair of 
points the largest (3 value, 0max, for which the ^-neighborhood of the pair is empty. 
The motives are twofold: algorithmic and applied. From an algorithmic point of 
view this is an interesting problem for which no efficient solution was known before. 
And from an applied perspective this has rich potential for the problem of character 
recognition since the spectra of a set of points can be construed as its signature. 
Moreover, any /^-skeleton can be extracted from the ^-spectrum in linear time. 

Our algorithms are in log n+K) time for pmax > 1 and in 0(n^ log n+K) 

time for Pmax > 0. A point that lies in the ^-neighborhood of an edge is a witness 
to this edge. The quantity K is & count of the number of times that the edges of 
the initial (Delaunay or complete) graph on P are witnessed for P = 00 . 
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The rest of the chapter is organized as follows. In the next section we present 
algorithms based on the techniques of range searching and sweepline. In the two 
following sections we discuss respectively how to compute efficiently the A:/?-spectrum 
and the additively weighted /5-spectrum. We conclude in the fifth and final section. 


4.2 /3-spectrum 


In this section we propose two different algorithms for computing the /3-spectrum 
of P. 

The main idea of both the algorithms is to start with a sufficiently large super 
graph on P and compute the witness set Wg for each edge e, setting /5 = oo. For 
/? > 1 we start with Delaunay triangulation on P since it is a super-graph of the 
Gabriel graph (corresponds to /? = 1). For /3 > 0, we start with the complete graph 
on P. We discuss only the first case as the second case can be handled similarly. 

We begin with the key observation that all the witnesses of an edge are contained 
in the lune that corresponds to 0 = oo. This an open infinite strip bounded by the 
lines normal to the edge at its end-points. We assume without loss of generality 
that no edge is horizontal, since all such edges can be handled separately in a 
preprocessing step in 0(nlog n) time. 

Let T be a transformation that maps a point p, with coordinates (a, b), into the 
line Tp with equation y = ~ax -f b and a line I with equation y = mx c into the 
point Ti with coordinates (— m, c) (see Fig. 24). This is known as a pont-line duality 
transform. 

Under this point-line duality transformation T, an infinite strip in the primal 
plane maps to a vertical line segment in the dual plane. So a point p G Wg if and 
only if Tp intersects Tg (see Fig. 25). The problem of computing We for an edge e, 
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T 

• ^ 

P 

Figure 24: Geometric dual transformation T 

therefore, reduces to computing all the lines that intersect a vertical line segment 
corresponding to e in the dual plane. 

We consider two different approaches to the problem. 



4.2.1 The Range-Searching Approach 

In this approach we attempt to preprocess the set of lines L such that given a vertical 
query segment q, we can report the lines of L that intersect q in 0(log n + We) time, 
where We is the number of intersecting lines. 

A naive approach would be to locate one end-point of the query line-segment in 
the arrangement A of the lines L, and move to the other end-point, recording all 
intersections with the the lines of the arrangement. This turns out to be inefficient 
because the complexity of an individual cell may preclude a constant-time transition 
to an adjacent cell. But this is precisely what we want, and can be achieved at the 
cost of some extra preprocessing. 

In order to reduce the complexity of an individual cell so that we can move from 
one cell to an adjacent one in constant time, we first trapezoidize the cells of the 
arrangement by drawing vertical lines through the vertices of the arrangement. Each 
such line is terminated by the first line in A that it meets, if any at all. This can 
be done using the sweepline technique in 0{n^ log n) time. The resulting geometric 
structure is shown in Fig. 26. Let us name this planar subdivision A'. 




Figure 26; New planar subdivision A' 


Though each cell now has the geometric appearance of a trapezium, there can 
be many vertices on its non-parallel sides (see Fig. 27 ) in which case it would still 
be very complex. Such a trapezium has an upper (resp. lower) anomaly if its upper 
(resp. lower) side consists of more than three vertices. 

For an edge e 6 . 4 , let /I = {ai,a2, .. be the set all vertical edges in A! 
that are incident on e from above and B = {f>i, ^2, ■ ■ • j be the set of those edges 
that are incident from below. 

For j in [1 ... i; - 1] let {oi, aj+i, . . . , %•} be the subset of A whose end-points on 
e lie strictly between those of bj and 6^4.1. We extend the edges Oj+i, Oj+s, . . . , Ci+i, 
where I < {j - i) (see Fig. 28 ) up to an edge of A lying immediately below e. While 
this removes the upper anomalies of all the trapezium that hang on e, it may create 
upper anomalies for the ones that lie immediately below. 

To handle all the regions uniformly, we process A' by levels of the original 
arrangement A from the topmost level to the bottommost one, removing all upper 
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anomalies. And again from the bottom to the top, removing all lower anomalies. 
The data structure underlying the resulting planar subdivision A!' is called a hive- 
graph H [Cha86]. 

The time and space complexities of this preprocessing step are in 0(n2 log n) 
and O(n^) respectively. 

We preprocess A" for point location in 0{n^ log n) time [EGS86]. We index the 
regions of A", and store these regions in an array R such that the array index of a 
region is identical to its region index. Along with a region we also store the index 
of the region just above it. Given a query segment g, we first locate its bottom 
end-point in A" in (9(log n) time; we then use its region index to locate it in the 
array R and determine the index of the region just above. We repeat this till we hit 
the top end-point of q. 

We claim Ahe following result. 


Lemma 4.1 A set of n lines can be preprocessed in 0(n^ log n) time using O(n^) 
space such that given a vertical query segment we can report all lines that intersect 
it in 0(log n + Wg) query time, where Wg is the size of the output. 


It follows from the above lemma that if we are given a set of n lines L and a 
set of 0(n) vertical line segments V, the intersections of the lines in L with the line 
segments in V can be computed in 0(n^ log n + K) amortized time, where JsT is a 
count of all such intersections. 

The complexity of this algorithm is prima facie worse than than that of the 
straight forward algorithm. But we can whittle away a y/n factor from the above 
complexity by preprocessing groups of lines of size y/n. The details are in the lemma 
below. 
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Lemma 4.2 Given a set L of n lines and a set V of 0(n) vertical segments, 
computing the intersections of the vertical segments and the lines in L can he done in 
iQg tijne using 0{n) space, where K is a count of all such intersections. 


Proof: We divide the problem into m subproblems, each consisting of n/m lines 

and all the vertical segments in P. Solve each subproblem by the algorithm above. 
Summing up the time complexities of the individual subproblems, we find that the 
optimal time complexity of log n + K) is achieved for m = The space 

required is 0{n) since we are considering yjn lines at a time. | 


Knowing the potential witness set of each edge in DT, we can compute the praax 
values for each edge e in time proportional to the size of W^. Therefore ^rnax value 
for all the edges of DT can be computed in time proportional to K. Thus we have 
the following result. 


Theorem 4.3 The ^-spectrum of a set of points P for p > 1 can be computed in 
0{n^l'^ log n + K) time using 0{n) space, where K is a count of the total number of 
times edges in the DT are witnessed for P = oo. 


Similarly, we can obtain the following result. 


Theorem 4.4 The P-spectrum of a set ofn points P for P>0 can be computed in 
0(n^log n) time. 


We note that in this case we have to query Oin'^) vertical line segments and 
hence dividing the lines into groups of size ^/n does not help. 
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4.2.2 Sweepline Approach 


We can avoid a complicated data structure like the hive-graph, and at the same time 
take advantage of the fact tliat we know all query segments beforehand by using a 
sweepline technique. 

In the sweepline status £ we store the lines in L according to a total ordering 
of their intersections. The event queue Q is an union of the set of vertical segments 
V and the intersections of the lines in L, sorted by their x-values. 

The algorithm works as follows. As we sweep the structure, if the next event 
corresponds to an intersection of two lines we update £ by reordering their positions 
in the intersection-order with the sweepline. If the event is a vertical line segment, 
we report all the lines in £ that intersect it. This can be done in 0(log n -l- We) 
time, where We is the number of lines that intersect the vertical line segment. 

As in the range searching technique we also have the following two theorems 
w'hich follow from the observations above. 

Theorem 4.5 The 0 -spectrum of a set of points P for /? > 1 can be computed in 
log n + K) time using 0{n) space, where K is a count of the total number of 
times the edges in DT are witnessed for /? = oo. 

Theorem 4.6 The 0-spectrum of a set of points P for 0 >0 can be computed in 
0{n^ log n-f- AT) time using 0(n®) space, where K is a count of the total number of 
times the edges of complete graph are witnessed for 0 = oo. 


4.3 /c/3"Spectrum 


As we noted in an earlier chapter, there is a natural generalization of the notion of 
a /^-skeleton to that of a /cj5-skeleton obtained by relaxing the emptiness criterion 
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of a ^-lune to let it contain at most k - 1 points. Similarly, the A;/3-spectrum is 
a generalization of the /?-spectrum. The /c, 5-spectrum of a set of n points P is a 
weighted geometric graph G. Each edge e e G assigned with largest ^ value /^max) 
for which the /^-neighborhood of the edge contain at most - 1 points. 

To compute the spectrum for /? > 1, we start with the A:-Delaunay graph {k~ 
DG) [SC90], known to be a super graph of a A:/?-skeleton. To compute the spectrum 
for ^ values greater than or equal to 0, we start with the complete graph. 

Similar to previous methods we compute the set of witnesses for each edge 
e e A:-DG. For each witness w G We compute the p value at which w becomes 
witness of the edge e. The P,nax of the edge e is the kth smallest of these P values, 
ii\We\ > k otherwise P^ax = oo. Therefore, we claim the following result. 


Theorem 4.7 The kp -spectrum of a set P of n points in the plane for P > I, can 
be constructed in 0{kn^l^ log n + K) time using 0{n) space, where K is a count of 
the total number of times the edges in k-DT are witnessed for /? = oo. 


Similarly one can establish the following theorem. 


Theorem 4.8 The kP -spectrum of a set P of n points in the plane for P>0, can 
be constructed in 0{rP log n + A') time using 0{rP) space, where K is a count of 
the total number of times the edges of the complete graph are witnessed for p = oo. 


4.4 Weighted /^-Spectrum 

Let Wi be the additive weight of the point pi G P. Each weighted point pi can 
be treated as a circle Bj, centered at pu with radius equal to its weight Wi. The 
weighted distance du,ipi,x), between a point pi G P and an arbitrary point x, is 
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d{pi,x) Wi if 2 ^ is outside £i,-, and zero otherwise. The weighted distance dyj{pi p ) 
between any two points € P is the distance between Bi and Bj. Let D denote 
the set of circles Bi . 

The additively weighted .5-spcctrum is a weighted geometric graph, each edge e 
being \\eighted with the laigest /3 value, Pman for which the /^-neighborhood of e is 
not intersected by any weighted circle. 

In the rest of this section we discuss an algorithm for computing the weighted 
/5-spectrum. 


4.4.1 An Algorithm 

To compute the spectrum for /5 > 1 we start with a super graph which is the 
weighted Delaunay graph (WDG) [Mir93]. For p > 0, we start with the complete 
graph. We discuss only the first case. 

For each edge e € WDG, we compute the set of witnesses for /? = oo. The 
previous algorithms do not generalize in a straight forward manner, because a circle 
in the primal plane does not have a dual that is a simple object. However, a simple 
observation helps us get around this problem. To explain this, we need the help of 
some notations. The symbol Sfrpj denotes the luneooipi,Pj)\ we will let denote 
the smallest infinite strip which encloses the weighted circle Bk and is parallel to Se- 

The following observation is very straightforward. 


Observation 4.1 A point p* is witness to an edge e for some value of P if and only 
if SI intersects S^. 


The standard point-line dual transformation T maps a non-vertical infinite strip 
into a vertical line segment. Moreover, two parallel non-vertical infinite strips 
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Figure 29: Illustration of Observation 1. 

intersect each other in the primal plane if and only if their corresponding vertical 
line segmentvS overlap in dual plane. In other words, a point pk is witness of an edge 
e for some value of 0 if and only if the vertical line segments and Ts^ overlap 
each other in the dual plane. We use the simpler notations and /g for and 
Ts^ re.spectively. In terrn.s of this simpler notation the above observation may be 
restated as below: 


Observation 4.2 A point Pk is witness to an edge e for some value of 0 if and only 
if intersects /«. 


Therefore, to compute the set of witnesses of an edge e € WDG, we find all the 
vertical line segments /| that intersect the vertical segment in the dual plane. The 
question now is how efficiently can we get all the vertical segments J|’s for each e G 
WDG ? The following observations provide an insight. Their proofs are immediate 
from the definition of the dual transform. 
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as follows; if ch i« the toi) end point of /f, delete If from since If is below 
El 2 - Otheiwise it is the bot.toni end point of If] insert If in Cei^i since EI 2 is in 
If. Siinilciily npdntc Cei^ ^'ith respective to the point There is no chcinge in the 
interval lists, Cgh and C/i/.,. 

When the sweeplinc £ coincides with 1^, we report all the intersecting intervals. 
This is done as follows. \Ye find the elementary interval El e C such that the top 
end point of is in BI. The interval /« intersects all the intervals in Cei- This 
takes care of all intervals, whose top end point is above EL But the top end-points 
of some intersecting intervals may start below that of El. These can be reported by 
going down from the top end-point to the bottom end-point of in the sweepline 
status £. During this, if the top end-point of any interval is encountered we report 
the interval corresponding to this top end-point. 

It is eavSy to determine the events arising out of the adjacency of a pair of points 
on the sweepline status £. This is the x-value corresponding to the slope of a 
common tangent to the weighted circles (in the primal plane) of the intervals (in 
the dual plane) to which these end-points belong. This can be done in 0(1) time 
for each event (see Fig. 32). 




Figure 32: Finding an event. 
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Lemma 4.9 Ihc xiritiie.ss snt for all the edges e 6 WDG can be computed inO(n'^ log nH- 
K) time using 0{n^) space, for ft = oo. 

Proof: At any time, the swc('{)line status C contains 2n end points of n intervals. 

C can be implemented simply by an array of size 2n. At x = 0, the sweepline status 
C can be initialized in 0{n log 7 ?.) time. Each pair of adjacent end points in C forms 
an elementary interval. For each elementary interval El e C we compute a set of 
intervals C/?/, and maintain tlu'in in a height-balanced tree. All this can be done in 
0(71^ log n) time. 

The initial events are determined by adjacent end-points in C. We compute 
these events by traversing the sweeplinc £, from one end to other in 0{n) time. We 
also add all tin; vmtical segment /« for all e €WDG. This initializes the event-queue; 
the complexity of this is in (){n log n). 

Let e be an event which is a change of order of q 2 and qz (see Fig. 31). The cost of 
locating these point in £, and hence the three consecutive elementary intervals Eli, 

Eh, and Elz such that [f/ 2 ,r/;j] = Rh is in O(log n). Updating Cgh can be done in 
0(log n) time, since we need to make at most two deletions and two insertions to 
the list Cei-i- I»^^oI'ting the next events of q\ and qz, qz and qz, and qz and qi into 
Q can be done in 0(log n) time. There are at most O(n^) events since an interval 
can intersect another interval at most constant number of times. Therefore the total 
cost of processing all such events is in 0(n^ log n) time. 

When the sweepline £ coincides with £, in 0(log n) time we locate the ele- 
mentary interval El E C such that the top end-point of £ is in El. We report all 
intervals in Cei and travel down the sweepline £ from the top end-point of £ to the 
bottom end-point of £, reporting intervals whose top end-point are encountered. 
This can be done in 0{we) time since there are at most 2we end points between the 
top and bottom end-points of £, where tn* number of intervals h intersected. The 
total cost of processing all the vertical segments £ for e €WDG, is in 0{n log n+K) 
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time, where K is a count of the number of times that the edges of the WDG are 
witnessed for /? = cxd. | 

The complexity of the above algorithm is prima facie worse than the brute-force 
algorithm which tests a circle for intersection with each of the 0(n) infinite strips, 
corresponding to each edge of the WDG! For the complexity of this is in C>(n), and 
therefore in Oin^) for n circles. However, if we look at the above algorithm very 
carefully we notice that the problem lies in considering all the intervals in one go. 
The cost of maintaining the sweepline status of n intervals is high and dominates the 
algorithm. So, similar to previous algorithms here also we balance the processing 
time for intervals and finding the witnesses to improve the time complexity. One 
can easily verify the following lemma. 


Lemma 4.10 The set of witnesses for each edge e e WDG can be computed in 
0{n^/'^ log n + K) time using 0{n) space, where K is is a count of the total number 
of times edges in WDG are witnessed for /5 = oo. 

This gives us the following theorem. 


Theorem 4.11 For /? > 1, the weighted ^-spectrum of P can be constructed in 
0{n^l'^ log n + K) time using 0{n) space, where K is is a count of the total number 
of times edges in WDG are witnessed for p = oo. 

Similarly, the weighted /3-spectrum of P, for /3 > 0, can be constructed by 
starting with the complete graph on P as the super graph and find set of witnesses 
for each edge. The following theorem is easily proved. 

Theorem 4.12 The weighted P -spectrum of a set P of n points in the plane, for 
/3 > 0, can be constructed in 0{nP log n -f K) time using 0{rP) space, where K is 
is a count of the total number of times edges are witnessed for P = oo. 
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4.5 Conclusion 


In this chapter we have proposed efficient algorithms for computing the /5-spectrum, 
A: /5-spectrum and, weighted /5-spectrum in the L 2 metric. It would be interesting to 
implement these algorithms and test them on some practical data. 



Chapter 5 


/3- Skeletons in Higher Dimensions 


5.1 Introduction 


In this chapter we extend the definition of a ^^-skeleton to higher-dimensions. We 
show that the size of a lune-based /5-skeleton for /? > 2 in higher dimensions is linear. 
We also present an efficient algorithm for computing a /5-skeleton for /5 > 2. 

The rest of the chapter is organized as follows. In the next section we define the 
two types of /5-skeletons - lune and sphere-based. In the following section, we prove 
an upper bound on the size of a /5-skeleton for P > 2. We next present an efficient 
algorithm for computing a /5-skeleton for p > 2. We conclude in the third and final 
section. 


5.2 Definitions 

A /5-skeleton is a geometric graph, obtained by joining pairs of points whose P- 
neighborhoods are empty. The neighborhood size depends on the parameter /5, and 
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are of two types. 

The sphere-based /^-neighborhood of a pair of points {pi,Pj) is defined as follows: 
for P 1, it is the union of all open spheres of radius P\pipj\/2, passing through pi 
and pj. For p € [0, 1], it is the interior of the intersection of all open spheres of 
radius |pipj|/(2/3), passing through pi and pj. 

The lune-based /^-neighborhood of a pair {jpupp of points is defined as follows; 
for p > I, it is the interior of the intersection of two spheres of radius P\'p^\/2, 
centered at the points (1 - P/2)pi -l- {P/2)pj and {P/2)pi -f (1 - PI2)pj, respectively. 
For P G [0, 1], the lune-based neighborhood is identical with the sphere-based 
neighborhood. 


5.3 Algorithm 

5.3.1 Main Ideas 

The essential technique is to compute a small super graph of a /3-skeleton and then 
prune the edges of this super graph that do not belong to the /3-skeleton. 

To accomplish the first step we use the powerful region approach of Yao [Yao82]. 
This gives us a super graph, called the Geographical Neighborhood Graph (GNG, 
for short). 

In the section below, we discuss an improved version of Yao’s approach. To aid 
the intuition we discuss first the two-dimensional case and then indicate how the 
result can be extended to higher dimensions. 
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5.3.2 The two-dimensional case 

The ke}’ idea of this approach is to decompose the neighborhood of a point pi into 
a finite number of sectors by drawing rays from Pi. Each sector is said to be a 
narrow region and is characterized by the fact that if pj and pk are two points in 
a sector then dist{pj,pk) < max{dist{pi,pj),dist{pi,pk)). This latter inequality can 
be guaranteed by letting the angle of each sector be < 60°. The nice consequence is 
that in each sector only the edges connecting pi to its nearest neighbors are edges 
of the super graph. 

Our observation is that the upper bound of 60° degrees is a function of the 
parameter value (5. And though the inequality characterizing a narrow region need 
not hold when we raise the upper bound, the following is true. 



Figure 33: Region approach 


Lemma 5.1 For a given point pi, let the angle of a sector R beO = arccos(l/ V^)- 
If more than one point in R is nearest to pi, then none of the edges, connecting pi 
to its nearest neighbors is in a ^-skeleton for P > 2. 

Proof: Let the points Pj and Pk be two nearest neighbors of the point pi in the 

region R (see Fig. 33). Let the angle LpjPiPk = The points C,- and Cj are the 
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centers of circles whose radius is /3|pJpj|/2 and pass through pi and pj respectively, 
as shown in Fig. 34. Elementary geometry shows that the point pk is inside the 
,d-neighborhood of pi and Pj\ for it is trivial that CjPj is longer than CjPk- And the 
fact that Pj € /3-neighborhood (pi,pjfc), follows from the upper bound on the angle a. 

I 


Therefore, we divide the region around a point pi by rays that emanate from pi, 
making angles of (m- 1) arccos(l/v/2^) with the rc-axis, where l<m<6if0<72 
and 1 < m < 5 otherwise. We number the regions between two successive rays 
(m - 1) arccos(l/\/5^) and m arccos(l/\/2^) by m and label the region bounded 
by these as i?m(Pi)- 

An important consequence of the above lemma is that we get a super graph of 
linear size for /3 > 2. 

To construct the GNG, for each point p, we have to find its nearest neighbors in 
each sector jR,-. We show how to do this for the region Ri (the remaining regions can 
be handled similarly). We rotate the y-axis by an angle 90° - 9. This is equivalent 
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to transforming each point p into a new point p' by the following transformation. 


Px' 


1 -cot(0) 


Px 

. Py' . 


_ 0 l/sin(0) 


.Py . 



Figure 35: Divide the points into cells. 

We sort the new points by their X-coordinates, and divide them into m con- 
secutive subsets of equal size. We sort the points in each subset by their new 
F-coordinates and divide each subset into m consecutive subsets of equal size as 
shown in Fig. 35. For each cell Cj, we compute the following four attributes: 

^max(ci) = max{px'\p e a} 

Ymaxici) = rnax{py>\p € q } 

XminiCi) - min{pa^\p € Ci} 

Xminici) = min{py'\p € Ci} 

With respect to p, the cells are divided into three classes. A cell is in class 
1, if all its points are in Riip)- If all its points are outside i?i(p), it is in class 2; 
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otherwise in class 3. To compute the nearest neighbors ofp in Ri{p) we use following 
procedure: cells of class 2 can be ignored. We compute the Voronoi diagram for each 
cell of class 1 in 0{n/m^ log{n/m^) time per cell, since the number of points in each 
cell is in 0{n/m^). We locate p in each cell’s voronoi diagram and find the nearest 
neighbor in each cell in O(log(n/m2)) time. We scan all the points in cells of class 
3 and find the nearest neighbor of p that is in Ri{p). 

Therefore, the total time required for computing nearest neighborhood in R^ 
for each point is in 0{n \og{n/m?) + nm^ \og{n/m?-) + ri^jm). Choosing m = 
(n/ logn)^/^, the time complexity is seen to be in log^/^ n). 


5.3.3 The higher-dimensional case 

To extend the above discussion to higher dimensions, we need the notion of a basis 
and a frame as introduced by Yao [Yao82]. We consider our metric space to be a 
d-dimensional real vector space. If B = is a basis of this vector space, 

then the convex cone, Conv{B) = {Eitri > 0 for all z} corresponds to what 
we called a sector in the two-dimensional case. A set of bases Bi, B 2 , ..., Bk form a 
frame, if U Cone{Bi) is the entire vector space. 

The angle of a cone is the maximum angle between any two vectors in this cone 
and the angle of the frame is the maximum angle taken over all the cones. 

Yao proved the interesting result that if ^ is any angle in the range 0 < V’ < tt, 
we can always construct a frame whose angle is less than ip. For our problem, we 
need to construct a frame whose angle is less than arccos(l/ V^). The independence 
of this angle with respect to dimension is a bit of a surprise. 

As in the 2-dimensional case, the following lemma identifies the possible edges 
of a /?-skeleton for P > 2. 
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Lemma 5.2 For a given point pi, let the angle of a cone R be bounded by 6 = 
arccos(l/v^)- If more than one point in R is nearest to pi, then none of the edges, 
connecting Pi to its nearest neighbors is in a p-skeleton for P > 2. 


Proof: Let the points pj and pk be two nearest neighbors of the point Pi in the 

region R. Let h be a plane passing through the points andp^t. The intersection 
of h and the ^-neighborhood(pi,Pj) is a 2-dimensional ^-neighborhood. Therefore, 
the point pk is in /?-neighborhood(pi,pj) from Lemma 5.1. Similarly, we can show 
that the point pj is in the /?-neighborhood(pi,pjfc). | 


As in the 2-diraensional case, we conclude that for fixed d, the size of the P- 
skeleton is in 0(n) for P >2. 

Given a frame T, the GNG can be computed using the method proposed by 
Yao in [Yao82]. Here is a brief sketch of his algorithm. Cnsider a basis Bi. We 
transform the coordinates of our points to this new basis Bi. Sort the points on the 
first coordinate and divide them into m consecutive subsets of equal size. In each 
suhset, sort the points on the second coordinate and divide these into m consecutive 
subsets of equal size. Repeat this processes for the remaining d—2 coordinates. This 
divides the points into subset of equal size in 0{dn log n) time. We preprocess 
each subset such that, given a query point p, we can find a nearest neighbor in 
0(log(n/m‘^)) time. For each point p, its nearest geographical neighbors can be 
computed as follows. Imagine the frame ^ to be placed at p. Consider a cone R. 
\\'ith respect to R, the mP cells are divided into three classes. A cell is in class 1, if 
all its points are in R. If all its points are outside R, it is in class 2; otherwise it is 
in class 3. To compute the nearest neighbors of p in R we use following procedure: • 
cells of class 2 can be ignored. We locate p in the Voronoi diagram of each cell in 
class 1 and find the nearest neighbor in 0{log{nfnnP)) time for this cell. We scan all 
the points in the cells of class 3 and find a nearest neighbor of p that is in R. This 
process, repeated for all the cones at p and all the input points, gives us the GNG. 
Following the analysis of [Yao82], we can prove the following result. 
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Lemma 5.3 Given a sat P of n points, a super graph of the ^-skeleton for p > 2 
can be computed in time, where a(d) = 

5.3.4 Pruning the GNG 

A brute-force method of thinning the super graph to the /3-skeleton is to check the 
status of each edge. This means checking if the lune corresponding to an edge is 
empty or not. Since this requires 0{n) time and there are 0(n) edges, the complexity 
is in 0(n^). 

To improve on this result, we consider all the spheres that determine the lunes of 
the edges in GNG. The resulting arrangement of spheres divides the space into 0(n‘^) 
cells. We can preprocess these hyperspheres, to obtain for each A:-face, 0 < k < d, 
the set of lunes that contain this A:- face. For each point p, locate the face / in 
the arrangement of hyperspheres and prune the all the lunes that contain /. That 
is, the problem of checking emptiness of lunes is reduced to point location in an 
arrangement of hyperspheres. 

Since an arrangement of hyperplanes is easier to work with, by an inversive 
transform we map the spheres to a set of hyperplanes in d + 1 dimensions. We 
briefly discuss the inversion transformation. 

An inversion is determined with respective to two parameters. These are the 
center and radius of inversion. Consider the planar case, with the center of inversion 
at the origin and radius of inversion Z > 0. A point p is mapped into another point p 
such that d{o,p) d{o,p') - L In other words, a vector in the direction 6 is mapped to 
another vector in the same direction but with inverted magnitude. Similarly, a circle 
which passes through the center of inversion is mapped to a line which does not pass 
through the center of inversion. Moreover, the interior of the circle is mapped to a 
half-plane as shown in the Fig. 36, 
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Figure 36: Inverse transformation in plane. 


The inversion of our hyperspheres to hyperplanes in a space that is one dimen- 
sional higher is carried out this way. Imagine our spheres to be embedded in a d- 
dimensional hyperplane in d+ 1-space. Let c be a center of inversion, lying outside 
this hyperplane. Consider a hypersphere, S', that passes through this center of 
inversion and a given sphere, S, that lies embedded in the d-dimensional hyperplane. 
The d-dimensional hyperplane, that S' inverts to with repect to c corresponds to S. 

Now the problem of point location in the arrangement of hyperspheres in d-space 
is reduced to the problem of point location in the arrangement of hyperplanes in 
(d -(- l)-space. Using the generalized binary search method for higher dimensions, 
proposed by Dobkin and Lipton [DL76], we can preprocess these hyperplanes for 
point location query. We briefly describe the method, first in 2-dimensions and then 
in higher dimensions. 

Given a set of n lines in the plane, these lines are preprocesed for point location. 
The projections of the intersections of these lines on the x-axis are sorted. These 
projected points divide the x-axis into O(n^) intervals, such that in each interval no 
two given lines intersect. Therefore, for each interval store the relative ordering of 
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the given lines. The time recgiiivd is in 0{n^ log n + n^n log n). The query can be 
aiKSwercKl in O(log n + log 7 ;’^) time. That is, 0(3 log n) time. 

We proceed similarly in {d+ l)-space. Project the intersc'ctions of all pairs of the 
giv<'n n hyperplamis on to ;i hypcuplane, which is distinct from given hyperplanes. 

In eacli region in the arrangemcuit of these O(n^) projected hyperplane, no two 
original hypcrplanes intersect. So for each region store the relative ordering of 
these hyperplanes. To locatr* the region in O(n^) hyperplanes in d-dimension, we 
recursively solve; the problem. 

Let P(7i, d +1) be the time required for preprocessing the n hypefplanes in 
(d + l)-dimension. We get the following recurrence relation for V{n,d): 

P(n, d + 1) = P(n^, d) + 0{ri^‘‘n log n) 

with V{n, 2) = 0{n^ n log n). 

L<!t Q{n, d + 1) the* time* r(;(iiiired for answering a query. We have the following 
r(‘curr(;ncc relation for Q{ii,d): 

Q(n, d + 1) = Q{n^, d) d- O(log n) 

with Q(n,2) = 0(log 71 ). 

The solutions to these r(;currences allow us to conclude the following. 

Lemma 5.4 A set ofn hyperplanes in {dA-\)-dimensio7i ca.n preprocessed in 0(n^ ^ h 

time so that each subsequent point location query can be answered in ^ + 

(d — 1)) log n) time. 


The complexity of th(; above algorithm is worse than the brute force algorithm 
for pruning the lunes. However, if we look at the above algorithm carefully we 
noticx; that the preprocessing t,ime for 0{n) hyperplanes is high and dominates the 
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algoiitlim, whciCES the (ju6iy time is low. Wo Ccin improve on the complexity by 
balancing the preprocessing and query times to derive the following result. 

Lemma 5.5 Given a GNG, the nonempty lunes can be eliminated in 
log n) time for P > 2, where a'{d) = (2^^+^ — 1)“^). 

Proof: We divide the problem into 0{n/m) subproblems, each subproblem con- 

sisting 0(m) hyperplanes and all the input points. For each of these subproblems, we 
preprocess the hyperplanes for point location in 0((n/m) log m) time. The 

total point location time is in 0(n^/m log m). Therefore the total time required to 
solve all the 0{n/m) subproblems is in 0[{n/m) log m-fn^/mlog m). The 

minimum time of 0{n^~°‘'^'^^ log n) is achieved for m = n°‘'^'^\ I 

The only remaining problem is to generate the lune list for each face. This 
can be done with the preprocessing time of hyperplanes [SC91c]. Hence follows the 
theorem. 

Theorem 5.6 A ^-skeleton for f3 > 2 can be computed in (log 

time. 


5.4 Conclusion 


In this chapter, we have generalized the notion of the ,0-skeleton to higher dimensions 
and presented an elSicient algorithm for computing a lune-based ^-skeleton for /3 > 2 
in any dimension. 


The problem of finding an optimal algorithm for computing ^-skeletons remains 
open. A recent interesting discovery is that in 2-dimension a circle-based /5-skeleton 
is subgraph of a Minimum Weight Triangulation on n points for /5 > 1.7317. It 
would be interesting to know if this result extends to higher dimensions. 



Chapter 6 


Conclusions 


In this thesis, we have attempted to further some of the of the algorithmic aspects 
of the class of geometric graphs known as /^-skeletons. The highlights of our contri- 
bution, in chapterwise order, are as follows. 

In the first chapter, we have have shown how to compute efficiently a /3-skeleton 
for any >0. The complexity of our algorithm is in log n) for /3 > 1 and in 

log n) for /3 in in the range [0, 1), for any Lp, where 1 < p < oo. In the Li 
and Loo metrics, the algorithms are in log n) for any /? > 0. 

We have also extended our algorithms to compute generalized /^-skeletons, namely, 
A;/?-skeleton and additively weighted /5-skeleton. The time complexity of computing 
the A:/3-skeleton is in 0{kn^^^ log n + k^n) for /? > 1 and in log n •+ n^k) for 

(3 in the range [0, 1). The time complexity of computing the weighted y3-skeleton is 
in log n) for /3 > 1 and in log n) for /? in the range [0, 1). 

In the metric Too, a /3-lune is simply a rectangle parallel to one of the axes. In 
the second chapter, we have exploited this property to present an optimal output- 
sensitive algorithm for computing a /3-skeleton for /3 > 2. The time required is in 
0{n log n -f- K), where K is size of the output graph. 
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In the third clicxptcr, \V6 htivc prnsented cin 6fRci6nt cilgorithrn. for computing 
the /3-spcctrum. Our algorithms are in 0{n^l'^ log n + K) for (3max > 1 and in 
0{ri^ log n + K) for Pmax ^ 0, where TC is a count of the number of times that the 
edges of the initial (Delaunay and complete respectively) graph are witnessed for 
/? = oo. 

We have also extended the notion of the /5-spectrum to A:/?-spectrum and addi- 
tively weighted /5-spectrum. We have described fast algorithms for computing these 
graphs. The time complexity for computing the /cyS-spectrum is in log n+K) 

for Prnax > 1 and in 0(n^ log n + K) time for pmax > 0, where K is a, count of the 
number of times that the edges of the initial (/c-Delaunay and complete respectively) 
graph are witnessed for /3 = oo. The algorithm for computing the weighted /?- 
spectrum takes 0(n^/^log n + K) time for Pmax > 1 and 0(n^ log n + K) time 
for Prnax > 0, where if is a count of the number of times that the edges of the 
initial (weighted Delaunay and complete respectively) graph on P are witnessed for 
f5 = oo. 

We have extended the notion of circle-based and lune-based /5-skeletons to higher 
dimensions. For fixed d, we have shown that the size of the lune based /5-skeleton 
in d dimension is linear for ^ > 2. We have also presented an efficient algorithm for 
computing a lune-based /5-skeleton for /5 > 2. 


6.1 Further Research Problems 

The research reported in this thesis opens up many new questions. 

A rather obvious one is to find more efficient algorithms for computing the 
various geometric graphs that we have discussed in this thesis, viz., /5-skeleton, A:/5- 
skeleton, weighted /5-skeleton etc. 

In the metrics Li and loo, computing the ;5-skeleton, for /5 in [1,2), in time lower 
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order than 0(n log n) is an interesting open question. It is also worth investigating 
if the output-sensitive algorithm presented in this thesis for computing the ^-skeleton 
can be generalized to higher dimensions for /? > 2. 

One really interesting avenue that we haven’t explored is the application of 
/9-skeletons to practical problems. We hope that somebody will take up this applca- 
tion challenge. Like, for example, using j0-skeletons for shape reconstruction and 
identification of real world objects. 

Recently, researchers have dicsovered an intriguing link between circle-based 
/9-skeletons and Minimum Weight Triangulations [Kei94, Yan95, CX96]. For a 
planar point set, a circle-based /9-skeleton is a subgraph of the Minimum Weight 
Triangulation for /9 > 1.17682. It would be worth investigating if this result 
generalises to higher dimensions. 

Another line of work that is tantalizingly open is the problem of recognising 
a /9-skeleton. Some preliminary work in this direction has been done by Bose et 
al [BLL96]. But a vast untapped field lies ahead. 



Appendix A 

Some Examples of /3-Skeletons 




p=1.00 


P=1.25 



p=1.50 
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